#include "Riostream.h"
#include <iostream>
#include <iomanip>
#include <string>

#include "TFile.h"
#include "TTree.h"
#include "TChain.h"
#include "TH1F.h"
#include "TH2F.h"
#include "TCanvas.h"
#include "TPad.h"
#include "TStyle.h"
#include "TROOT.h"
#include "TMath.h"
#include "TRandom.h"
#include "TFormula.h"
#include "TSystem.h"
#include "TEnv.h"

#include "CommandLineInterface.hh"
#include "Barrel.hh"
#include "Annular.hh"
#include "Germanium.hh"
#include "SRIMloader.hh"

using namespace TMath;
using namespace std;

fstream Nfile;

ClassImp(Barrel);
ClassImp(Annular);
ClassImp(Germanium);

SRIM srim;

// **********************************
// ********* PROGRAM START **********
// **********************************

int main(int argc, char* argv[]){

// -- -- -- Parameter definition -- -- --
	vector<char*> InputFiles;
	vector<char*> Histos;
	char* OutputFile = NULL;
	char* MBFile = NULL;
	char* cutfile = NULL;
	char* srimdir = NULL;
	char* SetFile = NULL;
	bool verbose = false;
	int tiltcorr = 0;
	int sim = 0;
	
	const Double_t pi = TMath::Pi();

	Double_t MBTheta[8][3][6];
	Double_t m[110];
	Double_t Ein, thetain, phiin, xin, yin;
	Double_t Pxin, Pyin, Pzin, Pin;
	Double_t ttarget, tBe;
	Double_t xpad[4], ypad[4], zpad[4];
	Int_t Apart, Atarget;
	Double_t PT, CDThreshold[4];	

	Double_t dEcor, Ecor;
	Double_t Edet[3], Efor[3][5], Q, EN, Eex[3][5];
	Double_t x, y, z;
	Double_t Px, Py, Pz;
	Double_t r, alpha, Rafter;
	Double_t theta, thetadet, thetap[3][5], thetain_coin, theta_cm[3][5];
	Double_t phi, phip[3][5], phiin_coin;
	Double_t Ein_coin, Ein_coin2;
	
	Float_t time[3][5];
	Double_t vcm, vout[2];
	Double_t beta[3][5], gamma[3][5];
	Double_t EGamDS;
	
	Int_t ID[3][5], par[3][5];
	Int_t l=0;
	Int_t Np, Nd, Nt;
	Int_t Np_FB, Nd_FB, Nt_FB;
	Int_t Np_CD, Nd_CD, Nt_CD;

// -- -- -- Some value assignments -- -- --

	// Masses
	m[0] = 1.00866491574*931.494e3;
	m[1] = 1.00782503207*931.494e3;
	m[2] = 2.01410177785*931.494e3;
	m[3] = 3.01604927767*931.494e3;
	m[4] = 4.00260325415*931.494e3;
	m[9] = 9.012182201*931.494e3;
	m[10] = 10.013533818*931.494e3;
	m[11] = 11.021657749*931.494e3;
	m[12] = 12.026920737*931.494e3; // 12Be
	m[20] = 19.999981315*931.494e3; // 20F
	m[21] = 20.999948951*931.494e3; // 21F
	m[22] = 22.002998815*931.494e3; // 22F
	m[23] = 22.994466904*931.494e3; // 23Ne
	m[107] = 106.905096820*931.494e3;

	// Input values, when runing the program
	CommandLineInterface* interface = new CommandLineInterface();

	interface->Add("-i", "inputfiles", &InputFiles);
	interface->Add("-o", "outputfile", &OutputFile);
	interface->Add("-M", "miniball angle", &MBFile);
	interface->Add("-C", "cut functions", &cutfile);
	interface->Add("-S", "SRIM files directory", &srimdir);
	interface->Add("-Set", "Setting file", &SetFile);
	interface->Add("-hi", "Histograms to be produced", &Histos);
	interface->Add("-t", "1 tiltcorr histos, 0 no tiltcorr", &tiltcorr);
	interface->Add("-v", "verbose", &verbose);
	interface->Add("-sim", "Set if you are analyzing a simulation",&sim);


	// input check
	interface->CheckFlags(argc, argv);

	if(InputFiles.size() == 0 || OutputFile == NULL){
		cerr<<"You have to provide at least one input file and the output file!"<<endl;
		exit(1);
	}
	if(Histos.size() == 0){
		cerr << "No histograms will be filled!"<<endl;
		exit(1);
	}
	if(srimdir == NULL){
		cerr << "SRIM directory is needed!"<<endl;
		exit(1);
	}

	cout<<"input file(s):"<<endl;
	for(unsigned int i=0; i<InputFiles.size(); i++){
		cout<<InputFiles[i]<<endl;
	}

	cout<<"output file: "<<OutputFile<< endl;

	// Loading settings
	TEnv* Settings = new TEnv(SetFile);
	

// -- -- -- Beam Parameters are set -- -- --
	// Energy
	Ein = Settings->GetValue("E.beam",2.85e3);
	Ein = Ein*11;
	// Angle
	thetain = Settings->GetValue("theta.beam",0.);
	thetain = thetain/180.*pi;
	phiin = Settings->GetValue("phi.beam",0.);
	phiin = phiin/180.*pi;
	// Offset
	xin = Settings->GetValue("x.beam",0);
	yin = Settings->GetValue("y.beam",0.);

// -- -- -- Target Parameters -- -- --
	ttarget = Settings->GetValue("t.target",10.);
	Atarget = Settings->GetValue("A.target",2);

	if (Atarget==12) m[12] = 12.000000000*931.494e3; // 12C

// -- -- -- Pad settings -- -- --
	for (Int_t ipad=0;ipad<4;ipad++){
		xpad[ipad] = Settings->GetValue(Form("x.pad.%d",ipad),0);
		ypad[ipad] = Settings->GetValue(Form("y.pad.%d",ipad),0);
		zpad[ipad] = Settings->GetValue(Form("z.pad.%d",ipad),0);
	}

// -- -- -- A for light particles, stopped in dE -- -- --
	Apart = Settings->GetValue("A.particle",2);

// -- -- -- Generating SRIM graphs -- -- --
	srim.Creategraphs(srimdir);

// -- -- -- Calculating beam energy and momentum at reaction -- -- --
	Rafter = srim.RT[11]->Eval(Ein);
	Ein = srim.ET[11]->Eval(Rafter - ttarget/(2.*cos(thetain)));
	
	Pxin = sqrt(2*m[11]*Ein)*sin(thetain)*cos(phiin);
	Pyin = sqrt(2*m[11]*Ein)*sin(thetain)*sin(phiin);
	Pzin = sqrt(2*m[11]*Ein)*cos(thetain);
	

// -- -- -- MINIBALL detector angles -- -- --
	TEnv *Angels = new TEnv(MBFile);
	for (int i=0; i<8; i++){
		for (int j=0; j<3; j++){
			for (int k=0; k<6; k++){
				MBTheta[i][j][k] = Angels->GetValue(Form("Theta.%d.%d.%d",i,j,k),0.0);
			}
		}
	}

// -- -- -- Loading all data -- -- --
	TChain* tr;
	tr = new TChain("caltr");
	for(unsigned int i=0; i<InputFiles.size(); i++){
		tr->Add(InputFiles[i]);
	}

	if(tr == NULL){
		cout << "could not find tree caltr in file " << endl;
		for(unsigned int i=0; i<InputFiles.size(); i++){
			cout<<InputFiles[i]<<endl;
		}
		return 3;
	}


// -- -- -- Parameters from the input files -- -- --
	vector<Barrel> *FBarrel = new vector<Barrel>;
	vector<Barrel> *BBarrel = new vector<Barrel>;
	vector<Annular> *FCD = new vector<Annular>;
	vector<Germanium> *Miniball = new vector<Germanium>;
	long long EbisTime;
	long long T1Time;
	long long SuperCycleTime;

	tr->SetBranchAddress("ForwardBarrel",&FBarrel);
	tr->SetBranchAddress("BackwardBarrel",&BBarrel);
	tr->SetBranchAddress("ForwardCD",&FCD);
	tr->SetBranchAddress("Miniball",&Miniball);
	tr->SetBranchAddress("EbisTime",&EbisTime);
	tr->SetBranchAddress("T1Time",&T1Time);
	tr->SetBranchAddress("SuperCycleTime",&SuperCycleTime);

// -- -- -- Output file is created -- -- --
	TFile* outfile = new TFile(OutputFile,"recreate");

	if(outfile->IsZombie()){
		return 4;
	}

// -- -- -- Low energy threshold for the CD detectors (keV) -- -- --
	CDThreshold[0] = 500.;
	CDThreshold[1] = 1300.;
	CDThreshold[2] = 500.;
	CDThreshold[3] = 800.;

// -- -- -- Generating cut functions for particle identification -- -- --
	TFormula* For_cut[4][16];
	TFormula* Back_cut[4];
	TFormula* CD_cut[6][4];
	TFormula* CD_PT_cut[6][4];
	TFormula* Be_cut[4];

	TEnv *Cutfunc = new TEnv(cutfile);

	// Forward Barrel
	for (Int_t istrip=0;istrip<16;istrip++){
		for (Int_t icut=0;icut<4;icut++) For_cut[icut][istrip] = new TFormula(Form("For_cut_%d_%d",icut,istrip),Cutfunc->GetValue(Form("For.%d.%d",icut,istrip),"0*x"));
	}
	// 11 Beryllium
	Be_cut[0] = new TFormula("Be_cut_0",Cutfunc->GetValue("For.11","0*x"));
	Be_cut[1] = new TFormula("Be_cut_1",Cutfunc->GetValue("For.12","0*x"));

	// Backward Barrel
	for (Int_t icut=0;icut<4;icut++) Back_cut[icut] = new TFormula(Form("Back_cut_%d",icut),Cutfunc->GetValue(Form("Back.%d",icut),"0*x"));

	// CD
	for (Int_t idet=0;idet<4;idet++){
		for (Int_t icut=0;icut<6;icut++){
			CD_cut[icut][idet] = new TFormula(Form("CD_cut_%d_%d",icut,idet),Cutfunc->GetValue(Form("CD.%d.%d",icut,idet),"0*x"));
			CD_PT_cut[icut][idet] = new TFormula(Form("CD_PT_cut_%d_%d",icut,idet),Cutfunc->GetValue(Form("CD_PT.%d.%d",icut,idet),"0*x"));
		}
	}
	// 11 Beryllium
	Be_cut[2] = new TFormula("Be_cut_2",Cutfunc->GetValue("CD.11","0*x"));
	Be_cut[3] = new TFormula("Be_cut_3",Cutfunc->GetValue("CD.12","0*x"));

// -- -- -- Prepare for vectorial histograms -- -- --
	// Forward Barrel
	
	// Backward Barrel
	
	// CD
	
	// Gamma

	
// -- -- -- Initiate histograms -- -- --

	// ** Single histograms **
	// Forward Barrel
	
	// Backward Barrel
	
	// CD
	
	// Gamma

	
	// ** Pad dependent histograms **
	for (Int_t i=0;i<4;i++){

		// ** Pad+Strip dependent histograms **
		for (Int_t j=0;j<16;j++){
	
		}

		// ** Pad+Angle dependent histograms **
		for (Int_t j =0;j<20;j++){
		}
	}

	// ** Angle dependent histograms **
	for (Int_t i=0;i<20;i++){
	
	}

	// ** CD dependent histograms **
	for(int i=0;i<4;i++){

		// ** CD+Strip/ring dependent histograms **
		for(int j=0;j<16;j++){
		}
	}

	// ** Gammas **
	// ** Cluster dependent **
	for (Int_t i=0;i<8;i++){

		// ** Cluster+Crystals dependent **
		for (Int_t j=0;j<3;j++){

			// ** Cluster+Crystal+Segments dependent**
			for (Int_t k=0;k<6;k++){
			}
		}
	}
// -- -- -- End initiating histograms -- -- --

// -- -- -- Start analysis -- -- -- 
	// Reset some values
	bool partfound;
	Double_t nentries = tr->GetEntries();
	Int_t nbytes = 0;
	Int_t status;

	Nd = 0;
	Np = 0;
	Nt = 0;
	Nd_FB = 0;
	Np_FB = 0;
	Nt_FB = 0;
	Nd_CD = 0;
	Np_CD = 0;
	Nt_CD = 0;
	tBe = 1;

// -- -- -- Start of loop over all data -- -- --
	for(int i=0; i<nentries;i++){
		// Reseting some parameters used for each event independently
		// (idet: 0=Forwardbarrel; 1=BackwardBarrel; 2=CD)
		// (ie: Event number in case of multihits)
		for (int idet=0;idet<3;idet++){
			for (int ie=0;ie<5;ie++){
				par[idet][ie] = 0;
				time[idet][ie] = 0;
				Eex[idet][ie] = 0;
				gamma[idet][ie] = 0;
				beta[idet][ie] = 0;
			}
		}

		// Checking event
		status = tr->GetEvent(i);
		if(status == -1){
			cerr<<"Error occured, couldn't read entry "<<i<<" from tree "<<tr->GetName()<<" in file "<<tr->GetFile()->GetName()<<endl;
			return 5;
		}
		else if(status == 0){
			cerr<<"Error occured, entry "<<i<<" in tree "<<tr->GetName()<<" in file "<<tr->GetFile()->GetName()<<" doesn't exist"<<endl;
			return 6;
		}
		nbytes += status;


		partfound =false;


// -- -- -- Analysing data from forward barrel -- -- --
		for(unsigned int j=0; j<FBarrel->size(); j++){
			// Checking the event
			if((*FBarrel)[j].GetMult() ==0) continue;

			if( (*FBarrel)[j].GetID() < 0 || (*FBarrel)[j].GetID() > 3 ){
				cerr<<"Error in entry "<<i<<": "<<j<<". forward barrel detector id is wrong: "<< (*FBarrel)[j].GetID()<<endl;
				continue;
			}

			if( (*FBarrel)[j].GetStripNr()[0] < 0 || (*FBarrel)[j].GetStripNr()[0] > 15 ){
				cerr<<"Error in entry "<<i<<": "<<(*FBarrel)[j].GetID()<<". forward barrel detector strip nr is wrong: "<<(*FBarrel)[j].GetStripNr()[0]<<endl;
				continue;
			}

			// Loop over all histo types
			for (Int_t ih=0;ih<Histos.size();ih++){
				partfound=true;

				// Time and energy of event stored
				time[0][j] = (*FBarrel)[j].GetTime();
				Edet[0] = (*FBarrel)[j].GetEdet();
				
				// Cylindrical coordinates
				x = 50*((*FBarrel)[j].GetStripPos()[0]-.5) + xin*cos(pi/2.*(*FBarrel)[j].GetID()) + yin*sin(pi/2.*(*FBarrel)[j].GetID());
				y = 29. + xin*sin(pi/2.*(*FBarrel)[j].GetID()) - yin*cos(pi/2.*(*FBarrel)[j].GetID());
				z = 8. + 3.125*((*FBarrel)[j].GetStripNr()[0] + gRandom->Rndm()); 
				// Adjusting coordinates if a detector is moved
				for (Int_t ipad=0;ipad<4;ipad++){
					if ((*FBarrel)[j].GetID()==ipad){
						x += xpad[ipad];
						y += ypad[ipad];
						z += zpad[ipad];
					}
				}

				r = sqrt(pow(y,2) + pow(x,2));
				alpha = atan(x/y);

				// Spherical coordinates
				theta = atan(r/z);
				phi = alpha + pi/2.*((*FBarrel)[j].GetID()+1);
				if (phi< 0) phi += 2*pi;
				if (phi> 2*pi) phi -= 2*pi;
				
				// Angle between particle and detector (norm)
				thetadet = acos(sin(theta)*cos(alpha));
				
				// Energy loss in carbon foil for alphas
				if (strcmp(Histos[ih],"alpha")==0){
					Rafter = srim.RDL[4]->Eval((*FBarrel)[j].GetRear());
					Efor[0][4] = srim.EDL[4]->Eval(Rafter + .2/cos(thetadet));
					Rafter = srim.RMylar[4]->Eval(Efor[0][j]);
					Efor[0][4] = srim.EMylar[4]->Eval(Rafter + 11.57/cos(thetadet));
				}
				
				// ** Filling histograms **
				
				// If silicon position is chosen
				if (strcmp(Histos[ih],"sipos")==0){
				}

				
// -- -- -- Events stopped in dE -- -- --
				// Check if particle is stopped
				PT = srim.EDL[Apart]->Eval(140./cos(thetadet));
				if (((*FBarrel)[j].GetEdet() < 400.) && ((*FBarrel)[j].GetRear() > 500.) && ((*FBarrel)[j].GetRear()<PT)) {
					
					// Values are set
					// (par: particle type: 1=proton; 2=deuteron; 3=triton; 4=Stopped particle (NPT))
					par[0][j] = 4;
					ID[0][j] = (*FBarrel)[j].GetID();
					thetap[0][j] = theta;
					phip[0][j] = phi;
					
					// ASSUMING ONE PARTICLE TYPE!!!!
					// Particle type given by Apart

					// Energy loss in carbon foil
					Rafter = srim.RMylar[Apart]->Eval((*FBarrel)[j].GetRear());
					Efor[0][j] = srim.EMylar[Apart]->Eval(Rafter + 11.57/cos(thetadet));
					
					// Energy loss in target
					if (Atarget==107){
						Rafter = srim.RT[20]->Eval(Efor[0][j]);
						Efor[0][j] = srim.ET[20]->Eval(Rafter + ttarget/(2.*cos(theta)));
					}
					else {
						Rafter = srim.RT[Apart]->Eval(Efor[0][j]);
						Efor[0][j] = srim.ET[Apart]->Eval(Rafter + ttarget/(2.*cos(theta)));
					}

					// Momentum vector of particle
					Px = sqrt(2*m[Apart]*Efor[0][j])*sin(theta)*cos(phi);
					Py = sqrt(2*m[Apart]*Efor[0][j])*sin(theta)*sin(phi);
					Pz = sqrt(2*m[Apart]*Efor[0][j])*cos(theta);

					// Energy of heavy fragment
					Q = m[11] + m[Atarget] - m[Apart] - m[11+Atarget-Apart];
					EN = (pow(Pxin-Px,2) + pow(Pyin-Py,2) + pow(Pzin-Pz,2))/(2*m[11+Atarget-Apart]);	
					Eex[0][j] = Ein + Q - Efor[0][j] - EN;
					
					// CM Angles
					vcm = (m[11]/(m[11]+m[Atarget]))*sqrt(2*Ein/m[11]);
					vout[0] = sqrt(2*Efor[0][j]/m[Apart])*cos(theta);
					vout[1] = sqrt(2*Efor[0][j]/m[Apart])*sin(theta);
					theta_cm[0][j] = atan(vout[1]/(vout[0]-vcm));
					if (theta_cm[0][j] < 0) theta_cm[0][j] += pi;

					// Relativistic factors
					gamma[0][j] = (Ein + (m[11] + m[Atarget] - m[11+Atarget-Apart] - m[Apart]) - Eex[0][j] - Efor[0][j] + m[11+Atarget-Apart])/m[11+Atarget-Apart];
					beta[0][j] = sqrt(1. - 1./(gamma[0][j]*gamma[0][j]));
					
					
					// ** Filling particle histograms **
					if (strcmp(Histos[ih],"part")==0){ 
					
						// Only detector 2 and 3 are used in forward direction !!!
						if (((*FBarrel)[j].GetID()==2)||((*FBarrel)[j].GetID()==3)){
							
						}
					}
				}
	
// -- -- -- Particle identification -- -- --
				if (((*FBarrel)[j].GetEdet() > 400.) && ((*FBarrel)[j].GetRear() > 500.)){
					Ecor = (*FBarrel)[j].GetEdet() + (1. - cos(thetadet))*(*FBarrel)[j].GetRear();
					dEcor = (*FBarrel)[j].GetRear()*cos(thetadet);
					
	
					// ** Protons **
					// (Analysed similar to stopped particles)
					if ((dEcor > For_cut[0][(*FBarrel)[j].GetStripNr()[0]]->Eval(Ecor)) && (dEcor < For_cut[1][(*FBarrel)[j].GetStripNr()[0]]->Eval(Ecor))){
						par[0][j] = 1;
						ID[0][j] = (*FBarrel)[j].GetID();
						Np++;
						Np_FB++;
						thetap[0][j] = theta;
						phip[0][j] = phi;

						Rafter = srim.RMylar[1]->Eval((*FBarrel)[j].GetEdet() + (*FBarrel)[j].GetRear());
						Efor[0][j] = srim.EMylar[1]->Eval(Rafter + 11.57/cos(thetadet));
						Rafter = srim.RT[1]->Eval(Efor[0][j]);
						Efor[0][j] = srim.ET[1]->Eval(Rafter + ttarget/(2.*cos(theta)));
						Px = sqrt(2*m[1]*Efor[0][j])*sin(theta)*cos(phi);
						Py = sqrt(2*m[1]*Efor[0][j])*sin(theta)*sin(phi);
						Pz = sqrt(2*m[1]*Efor[0][j])*cos(theta);
						Q = m[11] + m[Atarget] - m[1] - m[11+Atarget-1];
						EN = (pow(Pxin-Px,2) + pow(Pyin-Py,2) + pow(Pzin-Pz,2))/(2*m[11+Atarget-1]);
						Eex[0][j] = Ein + Q - Efor[0][j] - EN;
						
						vcm = (m[11]/(m[11]+m[Atarget]))*sqrt(2*Ein/m[11]);
						vout[0] = sqrt(2*Efor[0][j]/m[1])*cos(theta);
						vout[1] = sqrt(2*Efor[0][j]/m[1])*sin(theta);
						theta_cm[0][j] = atan(vout[1]/(vout[0]-vcm));
						if (theta_cm[0][j] < 0) theta_cm[0][j] += pi;
						
						// Relativistic factors
						gamma[0][j] = (Ein + (m[11] + m[Atarget] - m[11+Atarget-1] - m[1]) - Eex[0][j] - Efor[0][j] + m[11+Atarget-1])/m[11+Atarget-1];
						beta[0][j] = sqrt(1. - 1./(gamma[0][j]*gamma[0][j]));
				
					}
	
					// ** Deuterons **
					// (Analysed similar to stopped particles)
					if ((dEcor > For_cut[1][(*FBarrel)[j].GetStripNr()[0]]->Eval(Ecor)) && (dEcor < For_cut[2][(*FBarrel)[j].GetStripNr()[0]]->Eval(Ecor))){
						par[0][j] = 2;
						ID[0][j] = (*FBarrel)[j].GetID();
						Nd++;
						Nd_FB++;
						thetap[0][j] = theta;
						phip[0][j] = phi;
						
						Rafter = srim.RMylar[2]->Eval((*FBarrel)[j].GetEdet() + (*FBarrel)[j].GetRear());
						Efor[0][j] = srim.EMylar[2]->Eval(Rafter + 11.57/cos(thetadet));
						Rafter = srim.RT[2]->Eval(Efor[0][j]);
						Efor[0][j] = srim.ET[2]->Eval(Rafter + ttarget/(2.*cos(theta)));
						Px = sqrt(2*m[2]*Efor[0][j])*sin(theta)*cos(phi);
						Py = sqrt(2*m[2]*Efor[0][j])*sin(theta)*sin(phi);
						Pz = sqrt(2*m[2]*Efor[0][j])*cos(theta);
						Q = m[11] + m[Atarget] - m[2] - m[11+Atarget-2];
						EN = (pow(Pxin-Px,2) + pow(Pyin-Py,2) + pow(Pzin-Pz,2))/(2*m[11+Atarget-2]);
						Eex[0][j] = Ein + Q - Efor[0][j] - EN;
						
						vcm = (m[11]/(m[11]+m[Atarget]))*sqrt(2*Ein/m[11]);
						vout[0] = sqrt(2*Efor[0][j]/m[2])*cos(theta);
						vout[1] = sqrt(2*Efor[0][j]/m[2])*sin(theta);
						theta_cm[0][j] = atan(vout[1]/(vout[0]-vcm));
						if (theta_cm[0][j] < 0) theta_cm[0][j] += pi;
						
						gamma[0][j] = (Ein + (m[11] + m[Atarget] - m[11+Atarget-2] - m[2]) - Eex[0][j] - Efor[0][j] + m[11+Atarget-2])/m[11+Atarget-2];
						beta[0][j] = sqrt(1. - 1./(gamma[0][j]*gamma[0][j]));
					
					}
	
					// ** Tritons **
					// (Analysed similar to stopped particles)
					if ((dEcor > For_cut[2][(*FBarrel)[j].GetStripNr()[0]]->Eval(Ecor)) && (dEcor < For_cut[3][(*FBarrel)[j].GetStripNr()[0]]->Eval(Ecor))){
						par[0][j] = 3;
						ID[0][j] = (*FBarrel)[j].GetID();
						Nt++;
						Nt_FB++;
						thetap[0][j] = theta;
						phip[0][j] = phi;
						
						Rafter = srim.RMylar[3]->Eval((*FBarrel)[j].GetEdet() + (*FBarrel)[j].GetRear());
						Efor[0][j] = srim.EMylar[3]->Eval(Rafter + 11.57/cos(thetadet));
						Rafter = srim.RT[3]->Eval(Efor[0][j]);
						Efor[0][j] = srim.ET[3]->Eval(Rafter + ttarget/(2.*cos(theta)));
						Px = sqrt(2*m[3]*Efor[0][j])*sin(theta)*cos(phi);
						Py = sqrt(2*m[3]*Efor[0][j])*sin(theta)*sin(phi);
						Pz = sqrt(2*m[3]*Efor[0][j])*cos(theta);
						Q = m[11] + m[Atarget] - m[3] - m[11+Atarget-3];
						EN = (pow(Pxin-Px,2) + pow(Pyin-Py,2) + pow(Pzin-Pz,2))/(2*m[11+Atarget-3]);
						Eex[0][j] = Ein + Q - Efor[0][j] - EN;
						
						vcm = (m[11]/(m[11]+m[Atarget]))*sqrt(2*Ein/m[11]);
						vout[0] = sqrt(2*Efor[0][j]/m[3])*cos(theta);
						vout[1] = sqrt(2*Efor[0][j]/m[3])*sin(theta);
						theta_cm[0][j] = atan(vout[1]/(vout[0]-vcm));
						if (theta_cm[0][j] < 0) theta_cm[0][j] += pi;
						
						gamma[0][j] = (Ein + (m[11] + m[Atarget] - m[11+Atarget-3] - m[3]) - Eex[0][j] - Efor[0][j] + m[11+Atarget-3])/m[11+Atarget-3];
						beta[0][j] = sqrt(1. - 1./(gamma[0][j]*gamma[0][j]));
							
					}
				}
			}
		}
// -- -- -- End forward barrel -- -- --

// -- -- -- Backward barrel -- -- --
		for(unsigned int j=0; j<BBarrel->size(); j++){
			// Checking the event
			if((*BBarrel)[j].GetMult() ==0)
				continue;
			if( (*BBarrel)[j].GetID() < 0 || (*BBarrel)[j].GetID() > 3 ){
				cerr<<"Error in entry "<<i<<": "<<j<<". backward barrel detector id is wrong: "<<(*BBarrel)[j].GetID()<<endl;
				continue;
			}
			if( (*BBarrel)[j].GetStripNr()[0] < 0 || (*BBarrel)[j].GetStripNr()[0] > 15 ){
				cerr<<"Error in entry "<<i<<": "<<j<<". backward barrel detector strip nr is wrong: "<<(*BBarrel)[j].GetStripNr()[0]<<endl;
				continue;
			}

			// Loop over all histo types
			for (Int_t ih=0;ih<Histos.size();ih++){
				
				// Time and energy of event stored
				time[1][j] = (*BBarrel)[j].GetTime();
				Edet[1] = (*BBarrel)[j].GetEdet();
	
				// Cylindrical coordinates
				x = 50*((*BBarrel)[j].GetStripPos()[0]-.5) + xin*cos(pi/2.*(*BBarrel)[j].GetID()) + yin*sin(pi/2.*(*BBarrel)[j].GetID());
				y = 29. + xin*sin(pi/2.*(*BBarrel)[j].GetID()) - yin*cos(pi/2.*(*BBarrel)[j].GetID());
				z = 8. + 3.125*((*BBarrel)[j].GetStripNr()[0] + gRandom->Rndm()); 
				r = sqrt(pow(y,2) + pow(x,2));
				alpha = atan(x/y);
				
				// Spherical coordinates
				theta = atan(r/z);
				phi = alpha + pi/2.*((*BBarrel)[j].GetID()+1);
				if (phi< 0) phi += 2*pi;
				if (phi> 2*pi) phi -= 2*pi;
				
				// Detector angle
				thetadet = acos(sin(theta)*cos(alpha));
				
				
// -- -- -- Events stopped in dE -- -- --
				// Check if particle is stopped
				if (((*BBarrel)[j].GetEdet() < 200.) && ((*BBarrel)[j].GetRear() > 50.)){
					par[1][j] = 4;
					
					// ONLY PROTONS ANALYSED IN BACKWARD ANGLES!!!
					
					// Energy loss in carbon foil
					Rafter = srim.RT[1]->Eval((*BBarrel)[j].GetRear());
					Efor[1][j] = srim.ET[1]->Eval(Rafter + ttarget/(2.*cos(theta)));
					Px = sqrt(2*m[1]*Efor[1][j])*sin(pi-theta)*cos(phi);
					Py = sqrt(2*m[1]*Efor[1][j])*sin(pi-theta)*sin(phi);
					Pz = sqrt(2*m[1]*Efor[1][j])*cos(pi-theta);
					Q = m[11] + m[Atarget] - m[1] - m[11+Atarget-1];
					EN = (pow(Pxin-Px,2) + pow(Pyin-Py,2) + pow(Pzin-Pz,2))/(2*m[11+Atarget-1]);
					Eex[1][j] = Ein + Q - Efor[1][j] - EN;
					
					// CM Angles
					vcm = (m[11]/(m[11]+m[Atarget]))*sqrt(2*Ein/m[11]);
					vout[0] = sqrt(2*Efor[1][j]/m[1])*cos(pi-theta);
					vout[1] = sqrt(2*Efor[1][j]/m[1])*sin(pi-theta);
					theta_cm[1][j] = atan(vout[1]/(vout[0]-vcm));
					if (theta_cm[1][j] < 0) theta_cm[1][j] += pi;
					
					// Relativistic factors
					gamma[1][j] = (Ein + (m[11] + m[Atarget] - m[11+Atarget-1] - m[1]) - Eex[1][j] - Efor[1][j] + m[11+Atarget-1])/m[11+Atarget-1];
					beta[1][j] = sqrt(1. - 1./(gamma[1][j]*gamma[1][j]));
					
				}
	
				// Particle identification
				if (((*BBarrel)[j].GetEdet() > 200.)){
					Ecor = (*BBarrel)[j].GetEdet() + (1. - cos(thetadet))*(*BBarrel)[j].GetRear();
					dEcor = (*BBarrel)[j].GetRear()*cos(thetadet);
					
					// ** Protons **
					if ((dEcor > Back_cut[0]->Eval(Ecor)) && (dEcor < Back_cut[1]->Eval(Ecor))){
						par[1][j] = 1;

						// Energy loss
						Rafter = srim.RDL[1]->Eval((*BBarrel)[j].GetEdet() + (*BBarrel)[j].GetRear());
						Efor[1][j] = srim.EDL[1]->Eval(Rafter + .5/cos(thetadet));
						Rafter = srim.RT[1]->Eval(Efor[1][j]);
						Efor[1][j] = srim.ET[1]->Eval(Rafter + ttarget/(2.*cos(theta)));
						Px = sqrt(2*m[1]*Efor[1][j])*sin(pi-theta)*cos(phi);
						Py = sqrt(2*m[1]*Efor[1][j])*sin(pi-theta)*sin(phi);
						Pz = sqrt(2*m[1]*Efor[1][j])*cos(pi-theta);
						Q = m[11] + m[Atarget] - m[1] - m[11+Atarget-1];
						EN = (pow(Pxin-Px,2) + pow(Pyin-Py,2) + pow(Pzin-Pz,2))/(2*m[11+Atarget-1]);
						
						Eex[1][j] = Ein + Q - Efor[1][j] - EN;
						
						// CM Angles
						vcm = (m[11]/(m[11]+m[Atarget]))*sqrt(2*Ein/m[11]);
						vout[0] = sqrt(2*Efor[1][j]/m[1])*cos(theta);
						vout[1] = sqrt(2*Efor[1][j]/m[1])*sin(theta);
						theta_cm[1][j] = atan(vout[1]/(vout[0]-vcm));
						if (theta_cm[1][j] < 0) theta_cm[1][j] += pi;
						
						// Relativistic factors
						gamma[1][j] = (Ein + Q - Eex[1][j] - Efor[1][j] + m[11+Atarget-1])/m[11+Atarget-1];
						beta[1][j] = sqrt(1. - 1./(gamma[1][j]*gamma[1][j]));
	
					}
				}
			}
		}
// -- -- -- End backward barrel -- -- --

// -- -- -- Forward CD -- -- --
		for(unsigned int j=0; j<FCD->size(); j++){
			if( (*FCD)[j].GetID() < 0 || (*FCD)[j].GetID() >= 4 ){
				cerr<<"Error in entry "<<i<<": "<<j<<". backward barrel detector id is wrong: "<<(*FCD)[j].GetID()<<endl;
				continue;
			}

			time[2][j] = (*FCD)[j].GetTime();
			Edet[2] = (*FCD)[j].GetEdet();

			// ** Loop over all rings **
			for( unsigned int k=0;k<(*FCD)[j].GetRingNr().size();k++){
				if(verbose) cout << "ring nr " << (*FCD)[j].GetRingNr()[k] << " with energy " << (*FCD)[j].GetRingEnergy()[k] << endl;
				
				// Loop over histos types
				for (Int_t ih=0;ih<Histos.size();ih++){
					
					// ** Loop over all strips for each ring **
					for(unsigned int l=0;l<(*FCD)[j].GetStripNr().size();l++){
						if(verbose) cout << "strip nr " << (*FCD)[j].GetStripNr()[l] << " with energy " << (*FCD)[j].GetStripEnergy()[l] << endl;
	
						// ** Determining the angle of the particle
						// Values if reaction occurs in origin
						r = 9. + (gRandom->Rndm() + (*FCD)[j].GetRingNr()[k])*2.;	
						z = 63.;
						theta = atan(r/z);

						if (((*FCD)[j].GetStripNr()[l]>=4)&&((*FCD)[j].GetStripNr()[l]<12)){
							phi = .0593411*(4+ 2*((*FCD)[j].GetStripNr()[l]-4 + gRandom->Rndm())) - .71209;
						}
						else if ((*FCD)[j].GetStripNr()[l]<4){
							phi = .0593411*((*FCD)[j].GetStripNr()[l] + gRandom->Rndm()) - .71209;
						}
						else{
							phi = .0593411*((*FCD)[j].GetStripNr()[l]+8 + gRandom->Rndm()) - .71209;
						}
						
						// Moving reactionpoint to beamspot
						x = r*sin(phi) + xin*cos(pi/2.*(*FCD)[j].GetID()) + yin*sin(pi/2.*(*FCD)[j].GetID()); 
						y = r*cos(phi) + xin*sin(pi/2.*(*FCD)[j].GetID()) - yin*cos(pi/2.*(*FCD)[j].GetID()); 
						r = sqrt(pow(y,2)+pow(x,2));
						theta = atan(r/z);
						phi = atan(x/y);
						phi = phi + pi/2*((*FCD)[j].GetID()+1);
						if (phi>2*pi) phi -= 2*pi;
						if (phi<0) phi += 2*pi;
						
						
						// ** Check energy difference in ring and strip to identify true events **
						if (TMath::Abs((*FCD)[j].GetRingEnergy()[k] - (*FCD)[j].GetStripEnergy()[l]) < 500){
							
							// ** Particles stopped in dE **
							PT = srim.EDL[Apart]->Eval(500/cos(theta));
							if ((*FCD)[j].GetEdet() < CDThreshold[(*FCD)[j].GetID()]){
								par[2][j] = 4;
								thetap[2][j] = theta;
							
							// ASSUMING ONE PARTICLE TYPE!!!!
							// Particle type is given by Apart
							
							// Energy loss calculations
							if (Apart!=4){
								Rafter = srim.RDL[11+Atarget-Apart]->Eval((*FCD)[j].GetRingEnergy()[k]);
								Efor[2][j] = srim.EDL[11+Atarget-Apart]->Eval(Rafter + .7/cos(theta));
								if (Atarget==107){
									Rafter = srim.RT[20]->Eval(Efor[2][j]);
									Efor[2][j] = srim.ET[20]->Eval(Rafter + ttarget/(2.*cos(theta)));
								}
								else {
									Rafter = srim.RT[11+Atarget-Apart]->Eval(Efor[2][j]);
									Efor[2][j] = srim.ET[11+Atarget-Apart]->Eval(Rafter + ttarget/(2.*cos(theta)));
								}
							}
							
							// Momentum vector of light particle
							Px = sqrt(2*m[Apart]*Efor[2][j])*sin(theta)*cos(phi);
							Py = sqrt(2*m[Apart]*Efor[2][j])*sin(theta)*sin(phi);
							Pz = sqrt(2*m[Apart]*Efor[2][j])*cos(theta);
							
							// Energy of heavy fragment
							Q = m[11] + m[Atarget] - m[Apart] - m[11+Atarget-Apart];
							EN = (pow(Pxin-Px,2) + pow(Pyin-Py,2) + pow(Pzin-Pz,2))/(2*m[11+Atarget-Apart]);
							Eex[2][j] = Ein + Q - Efor[2][j] - EN;
							
							// CM Angles
							vcm = (m[11]/(m[11]+m[Atarget]))*sqrt(2*Ein/m[11]);
							vout[0] = sqrt(2*Efor[2][j]/m[Apart])*cos(theta);
							vout[1] = sqrt(2*Efor[2][j]/m[Apart])*sin(theta);
							theta_cm[2][j] = atan(vout[1]/(vout[0]-vcm));
							if (theta_cm[2][j] < 0) theta_cm[2][j] += pi;
      
							// Relativistic factors
							gamma[2][j] = (Ein + (m[11] + m[Atarget] - m[11+Atarget-Apart] - m[Apart]) - Eex[2][j] - Efor[2][j] + m[11+Atarget-Apart])/m[11+Atarget-Apart];
							beta[2][j] = sqrt(1. - 1./(gamma[2][j]*gamma[2][j]));
							
							
							// ** Coincidence events **
							
							// Loop over all barrel events
							for (Int_t ie=0;ie<5;ie++){
							
								// ** Coincidence with protons **
								if (par[0][ie]==1){
									Rafter = srim.RT[12]->Eval((*FCD)[j].GetStripEnergy()[l]);
									Efor[2][j] = srim.ET[12]->Eval(Rafter + ttarget/(2.*cos(theta)));
									Eex[2][j] = (pow(sqrt(2*m[12]*Efor[2][j])*cos(theta) + sqrt(2*m[1]*Efor[0][ie])*cos(thetap[0][ie]),2) + pow(sqrt(2*m[12]*Efor[2][j])*sin(theta) - sqrt(2*m[1]*Efor[0][ie])*sin(thetap[0][ie]),2))/(2*m[11]) - Efor[2][j] - Efor[0][ie] + (m[11] + m[2] - m[12] - m[1]);
		
									thetain_coin = atan((sqrt(2*m[12]*Efor[2][j])*sin(theta)*cos(phi) + sqrt(2*m[1]*Efor[0][ie])*sin(thetap[0][ie])*cos(phip[0][ie]))/(sqrt(2*m[12]*Efor[2][j])*cos(theta) + sqrt(2*m[1]*Efor[0][ie])*cos(thetap[0][ie])));
								}
									
								// ** Coincidence with deuterons **
								if (par[0][ie]==2){
									Rafter = srim.RT[11]->Eval((*FCD)[j].GetStripEnergy()[l]);
									Efor[2][j] = srim.ET[11]->Eval(Rafter + ttarget/(2.*cos(theta)));
									Eex[2][j] = (pow(sqrt(2*m[11]*Efor[2][j])*cos(theta) + sqrt(2*m[2]*Efor[0][ie])*cos(thetap[0][ie]),2) + pow(sqrt(2*m[11]*Efor[2][j])*sin(theta) - sqrt(2*m[2]*Efor[0][ie])*sin(thetap[0][ie]),2))/(2*m[11]) - Efor[2][j] - Efor[0][ie] + (m[11] + m[2] - m[11] - m[2]);
									
									thetain_coin = atan((sqrt(2*m[11]*Efor[2][j])*sin(theta)*cos(phi) + sqrt(2*m[2]*Efor[0][ie])*sin(thetap[0][ie])*cos(phip[0][ie]))/(sqrt(2*m[11]*Efor[2][j])*cos(theta) + sqrt(2*m[2]*Efor[0][ie])*cos(thetap[0][ie])));
									}
								
									// ** Coincidence with tritons
									if (par[0][ie]==3){
										Rafter = srim.RDL[10]->Eval((*FCD)[j].GetStripEnergy()[l]);
										Efor[2][j] = srim.EDL[10]->Eval(Rafter + .7/cos(theta));
										Rafter = srim.RT[10]->Eval(Efor[2][j]);
										Efor[2][j] = srim.ET[10]->Eval(Rafter + ttarget/(2.*cos(theta)));
										Px = sqrt(2*m[10]*Efor[2][j])*sin(theta)*cos(phi)+sqrt(2*m[3]*Efor[0][ie])*sin(thetap[0][ie])*cos(phip[0][ie]);
										Py = sqrt(2*m[10]*Efor[2][j])*sin(theta)*sin(phi)+sqrt(2*m[3]*Efor[0][ie])*sin(thetap[0][ie])*sin(phip[0][ie]);
										Pz = sqrt(2*m[10]*Efor[2][j])*cos(theta)+sqrt(2*m[3]*Efor[0][ie])*cos(thetap[0][ie]);
										Eex[2][j] = (pow(Px,2)+pow(Py,2)+pow(Pz,2))/(2*m[11]) - Efor[2][j] - Efor[0][ie] + (m[11] + m[2] - m[10] - m[3]);
										Ein_coin = (pow(Px,2)+pow(Py,2)+pow(Pz,2))/(2*m[11]);
											
											
										Ein_coin2 = Efor[2][j] + Efor[0][ie] - (m[11]+m[2]-m[10]-m[3]);
											
										thetain_coin = atan(sqrt(pow(Px,2)+pow(Py,2))/Pz);
										if (Px>=0&&Py>=0) phiin_coin = atan(Py/Px);
										if (Px>=0&&Py<0) phiin_coin = 2*pi + atan(Py/Px);
										if (Px<0&&Py>=0) phiin_coin = pi + atan(Py/Px);
										if (Px<0&&Py<0) phiin_coin = pi + atan(Py/Px);
											
									}
										
									// ** Coincidence with particles stopped in dE **
									if (par[0][ie]==4){
										Eex[2][j] = -Efor[2][j] - Efor[0][ie] + (m[11]+m[Atarget]-m[11+Atarget-Apart]-m[Apart]) + Ein;
											
										if (Apart!=4){
											Px = Pxin - (sqrt(2*m[Atarget+11-Apart-1]*Efor[2][j])*sin(theta)*cos(phi)+sqrt(2*m[Apart]*Efor[0][ie])*sin(thetap[0][ie])*cos(phip[0][ie]));
											Py = Pyin - (sqrt(2*m[Atarget+11-Apart-1]*Efor[2][j])*sin(theta)*sin(phi)+sqrt(2*m[Apart]*Efor[0][ie])*sin(thetap[0][ie])*sin(phip[0][ie]));
											Pz = Pzin - (sqrt(2*m[Atarget+11-Apart-1]*Efor[2][j])*cos(theta)+sqrt(2*m[Apart]*Efor[0][ie])*cos(thetap[0][ie]));
										
											Eex[2][j] = Ein - (pow(Px,2)+pow(Py,2)+pow(Pz,2))/(2*m[0]) - Efor[2][j] - Efor[0][ie] + (m[11] + m[Atarget] - m[11+Atarget-1-Apart] - m[Apart] - m[0]);
										}
									}
								}
							}
	
// -- -- -- Particle identification -- -- --
							if ((*FCD)[j].GetEdet() > CDThreshold[(*FCD)[j].GetID()]){
								Ecor = (*FCD)[j].GetEdet() + (1. - cos(theta))*(*FCD)[j].GetStripEnergy()[l];
								dEcor = (*FCD)[j].GetStripEnergy()[l]*cos(theta);
								
								// ** Protons **
								if ((dEcor > CD_cut[0][(*FCD)[j].GetID()]->Eval(Ecor)) && (dEcor < CD_cut[1][(*FCD)[j].GetID()]->Eval(Ecor))){
									par[2][j] = 1;
									Np++;
									Np_CD++;
	
									Rafter = srim.RT[1]->Eval((*FCD)[j].GetEdet() + (*FCD)[j].GetStripEnergy()[l]);
									Efor[2][j] = srim.ET[1]->Eval(Rafter + ttarget/(2.*cos(theta)));
									Px = sqrt(2*m[1]*Efor[2][j])*sin(theta)*cos(phi);
									Py = sqrt(2*m[1]*Efor[2][j])*sin(theta)*sin(phi);
									Pz = sqrt(2*m[1]*Efor[2][j])*cos(theta);
									Q = m[11] + m[Atarget] - m[1] - m[11+Atarget-1];
									EN = (pow(Pxin-Px,2) + pow(Pyin-Py,2) + pow(Pzin-Pz,2))/(2*m[11+Atarget-1]);
						
									Eex[2][j] = Ein + Q - Efor[2][j] - EN;
	
									// CM Angles
									vcm = (m[11]/(m[11]+m[Atarget]))*sqrt(2*Ein/m[11]);
									vout[0] = sqrt(2*Efor[2][j]/m[1])*cos(theta);
									vout[1] = sqrt(2*Efor[2][j]/m[1])*sin(theta);
									theta_cm[2][j] = atan(vout[1]/(vout[0]-vcm));
									if (theta_cm[2][j] < 0) theta_cm[2][j] += pi;
									
									gamma[2][j] = (Ein + (m[11] + m[Atarget] - m[11+Atarget-1] - m[1]) - Eex[0][j] - Efor[2][j] + m[11+Atarget-1])/m[11+Atarget-1];
									beta[2][j] = sqrt(1. - 1./(gamma[2][j]*gamma[2][j]));
									
								}
								
								// ** Deuterons **
								if ((dEcor > CD_cut[2][(*FCD)[j].GetID()]->Eval(Ecor)) && (dEcor < CD_cut[3][(*FCD)[j].GetID()]->Eval(Ecor))){
									par[2][j] = 2;
									Nd++;
									Nd_CD++;
									
									Rafter = srim.RT[1]->Eval((*FCD)[j].GetEdet() + (*FCD)[j].GetStripEnergy()[l]);
									Efor[2][j] = srim.ET[1]->Eval(Rafter + ttarget/(2.*cos(theta)));
									Px = sqrt(2*m[2]*Efor[2][j])*sin(theta)*cos(phi);
									Py = sqrt(2*m[2]*Efor[2][j])*sin(theta)*sin(phi);
									Pz = sqrt(2*m[2]*Efor[2][j])*cos(theta);
									Q = m[11] + m[Atarget] - m[2] - m[11+Atarget-2];
									EN = (pow(Pxin-Px,2) + pow(Pyin-Py,2) + pow(Pzin-Pz,2))/(2*m[11+Atarget-2]);
						
									Eex[2][j] = Ein + Q - Efor[2][j] - EN;
	
									// CM Angles
									vcm = (m[11]/(m[11]+m[Atarget]))*sqrt(2*Ein/m[11]);
									vout[0] = sqrt(2*Efor[2][j]/m[2])*cos(theta);
									vout[1] = sqrt(2*Efor[2][j]/m[2])*sin(theta);
									theta_cm[2][j] = atan(vout[1]/(vout[0]-vcm));
									if (theta_cm[2][j] < 0) theta_cm[2][j] += pi;
									
									gamma[2][j] = (Ein + (m[11] + m[Atarget] - m[11+Atarget-2] - m[2]) - Eex[0][j] - Efor[2][j] + m[11+Atarget-2])/m[11+Atarget-2];
									beta[2][j] = sqrt(1. - 1./(gamma[2][j]*gamma[2][j]));
									
								}
								
								// ** Tritons **
								if ((dEcor > CD_cut[4][(*FCD)[j].GetID()]->Eval(Ecor)) && (dEcor < CD_cut[5][(*FCD)[j].GetID()]->Eval(Ecor))){
									par[2][j] = 3;
									Nt++;
									Nt_CD++;
									
									Rafter = srim.RT[1]->Eval((*FCD)[j].GetEdet() + (*FCD)[j].GetStripEnergy()[l]);
									Efor[2][j] = srim.ET[1]->Eval(Rafter + ttarget/(2.*cos(theta)));
									Px = sqrt(2*m[3]*Efor[2][j])*sin(theta)*cos(phi);
									Py = sqrt(2*m[3]*Efor[2][j])*sin(theta)*sin(phi);
									Pz = sqrt(2*m[3]*Efor[2][j])*cos(theta);
									Q = m[11] + m[Atarget] - m[3] - m[11+Atarget-3];
									EN = (pow(Pxin-Px,2) + pow(Pyin-Py,2) + pow(Pzin-Pz,2))/(2*m[11+Atarget-3]);
						
									Eex[2][j] = Ein + Q - Efor[2][j] - EN;
	
									// CM Angles
									vcm = (m[11]/(m[11]+m[Atarget]))*sqrt(2*Ein/m[11]);
									vout[0] = sqrt(2*Efor[2][j]/m[3])*cos(theta);
									vout[1] = sqrt(2*Efor[2][j]/m[3])*sin(theta);
									theta_cm[2][j] = atan(vout[1]/(vout[0]-vcm));
									if (theta_cm[2][j] < 0) theta_cm[2][j] += pi;
									
									gamma[2][j] = (Ein + (m[11] + m[Atarget] - m[11+Atarget-3] - m[3]) - Eex[0][j] - Efor[2][j] + m[11+Atarget-3])/m[11+Atarget-3];
									beta[2][j] = sqrt(1. - 1./(gamma[2][j]*gamma[2][j]));
									
								}
							}
						}
					}
				}
			}
		}
// -- -- -- End CD -- -- --


// -- -- -- Gammas -- -- --
		// ** Loop over all MINIBALL events
		for(unsigned int j=0; j<Miniball->size(); j++){
			if((*Miniball)[j].GetCluID()<0 || (*Miniball)[j].GetCluID()>7){
				cout << "Error in entry "<<i<<": "<<j<<". Germanium detector Cluster id is wrong: "<<(*Miniball)[j].GetCluID()<<endl;
				continue;
			}

			// ** Loop over all crystals **
			for(unsigned int k=0; k<(*Miniball)[j].GetCrystal().size();k++){
				if( (*Miniball)[j].GetCrystal()[k].GetCryID() <0 || (*Miniball)[j].GetCrystal()[k].GetCryID() >2){
					cout << "Error in entry "<<i<<": "<<j<<". Germanium detector Crystal id is wrong: "<<(*Miniball)[j].GetCrystal()[k].GetCryID()<<endl;
					continue;
				}
				
				// ** Loop over all histos types **
				for (Int_t ih=0;ih<Histos.size();ih++){
	
	
					// Gammas in coincidence with particles
					for (Int_t idet=0;idet<3;idet++){
						for (Int_t ie=0;ie<5;ie++){
						
							// ** Gammas in coincidence with protons **
							if ((par[idet][ie]==1)&&(ID[idet][ie]!=0)){

								// Doppler correction
								EGamDS = gamma[idet][ie]*(*Miniball)[j].GetCrystal()[k].GetCore()*(1. - beta[idet][ie]*cos(MBTheta[(*Miniball)[j].GetCluID()][(*Miniball)[j].GetCrystal()[k].GetCryID()][(*Miniball)[j].GetCrystal()[k].GetSegID()[l]]));
							}

							// ** Gammas in coincidence with deuterons **
							else if ((par[idet][ie]==2)&&(ID[idet][ie]!=0)){
							
								// Doppler correction
								EGamDS = gamma[idet][ie]*(*Miniball)[j].GetCrystal()[k].GetCore()*(1. - beta[idet][ie]*cos(MBTheta[(*Miniball)[j].GetCluID()][(*Miniball)[j].GetCrystal()[k].GetCryID()][(*Miniball)[j].GetCrystal()[k].GetSegID()[l]]));
							}
									
							// ** Gammas in coincidence with Tritons **
							else if ((par[idet][ie]==3)&&(ID[idet][ie]!=0)){

								// Doppler correction
								EGamDS = gamma[idet][ie]*(*Miniball)[j].GetCrystal()[k].GetCore()*(1. - beta[idet][ie]*cos(MBTheta[(*Miniball)[j].GetCluID()][(*Miniball)[j].GetCrystal()[k].GetCryID()][(*Miniball)[j].GetCrystal()[k].GetSegID()[l]]));
							}
								
							// ** Gammas in coincidence with particles stopped in dE **
							else if ((par[idet][ie]==4)&&(ID[idet][ie]!=0)){
								// Doppler correction
								EGamDS = gamma[idet][ie]*(*Miniball)[j].GetCrystal()[k].GetCore()*(1. - beta[idet][ie]*cos(MBTheta[(*Miniball)[j].GetCluID()][(*Miniball)[j].GetCrystal()[k].GetCryID()][(*Miniball)[j].GetCrystal()[k].GetSegID()[l]]));
							
								// In coincidence with particles in backward barrel
								if ((idet==1) && ((ID[idet][ie]==2) || (ID[idet][ie]==3))){
							
								}
								else if (idet==0){
								}
							}
						}
					}
				}
			}
		}
// -- -- -- End gammas -- -- --

		
// -- -- -- Counter to show how far the program is -- -- --
		if(i%1000 == 0){
			cout<<setw(5)<<setiosflags(ios::fixed)<<setprecision(1)<<(100.*i)/nentries<<" % done\r"<<flush;
		}
	}

// -- -- -- Writing histograms -- -- --

	// Single histograms
	for (Int_t ih=0;ih<Histos.size();ih++){
	
		// Pad dependent histograms
		for (Int_t i=0;i<4;i++){
			
			// Pad+Strip dependent histograms
			for (Int_t j=0;j<16;j++){
			}
		}
		
		// Cd dependent histograms	
		for (Int_t i=0;i<4;i++){
				
			// CD+Strip/ring dependent histograms
			for (Int_t j=0;j<16;j++){
			}
		}
	
		// Cluster dependent histograms
		for (Int_t i=0;i<8;i++){
		
			// Crystals+Cluster dependent
			for (Int_t j=0;j<3;j++){

				// Segments+Crystal+Cluster dependent
				for (Int_t k=0;k<6;k++){
				}
			}
		}
	}
	
// -- -- -- Writing to outputfile -- -- --
	outfile->Close();
	delete tr;



// -- -- -- Number of detected particles pr. minute/second -- -- --
	printf("Total: Np = %d, Nd = %d, Nt = %d, t = %f min \n",Np/Histos.size(),Nd/Histos.size(),Nt/Histos.size(),EbisTime*25.e-9/60.);
	printf("Pr. Minute: Np = %f, Nd = %f, Nt = %f \n",Np/Histos.size()/(EbisTime*25.e-9/60.),Nd/Histos.size()/(EbisTime*25.e-9/60.),Nt/Histos.size()/(EbisTime*25.e-9/60.));
	printf("Pr. Minute in FB: Np = %f, Nd = %f, Nt = %f \n",Np_FB/Histos.size()/(EbisTime*25.e-9/60.),Nd_FB/Histos.size()/(EbisTime*25.e-9/60.),Nt_FB/Histos.size()/(EbisTime*25.e-9/60.));
	printf("Pr. Minute in CD: Np = %f, Nd = %f, Nt = %f \n",Np_CD/Histos.size()/(EbisTime*25.e-9/60.),Nd_CD/Histos.size()/(EbisTime*25.e-9/60.),Nt_CD/Histos.size()/(EbisTime*25.e-9/60.));


	Nfile.open("Npar.txt", fstream::out | fstream::app);

	Nfile << Np/(EbisTime*25.e-9/60.) << "	" << Nd/(EbisTime*25.e-9/60.) << "	" << Nt/(EbisTime*25.e-9/60.) << endl;

	Nfile.close();
	
	Nfile.open("Npar_FB.txt", fstream::out | fstream::app);

	Nfile << Np_FB/(EbisTime*25.e-9) << "	" << Nd_FB/(EbisTime*25.e-9) << "	" << Nt_FB/(EbisTime*25.e-9) << endl;

	Nfile.close();
	
	Nfile.open("Npar_CD.txt", fstream::out | fstream::app);

	Nfile << Np_CD/(EbisTime*25.e-9/60.) << "	" << Nd_CD/(EbisTime*25.e-9/60.) << "	" << Nt_CD/(EbisTime*25.e-9/60.) << endl;

	Nfile.close();
	
	Nfile.open("Time.txt", fstream::out | fstream::app);

	Nfile << EbisTime*25.e-9/60. << endl;

	Nfile.close();
	
	return 0;

}

